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Here we present an efficient and numerically stable procedure for compressing a correlation matrix 
into a set of local unitary single-particle gates, which leads to a very efficient way of forming the 
matrix product state (MPS) approximation of a pure fermionic Gaussian state, such as the ground 
state of a quadratic Hamiltonian. The procedure involves successively diagonalizing subblocks of the 
correlation matrix to isolate local states which are purely occupied or unoccupied. A small number of 
nearest neighbor unitary gates isolates each local state. The MPS of this state is formed by applying 
the many-body version of these gates to a product state. We treat the simple case of compressing 
the correlation matrix of spinless free fermions with definite particle number in detail, though the 
procedure is easily extended to fermions with spin and more general BCS states (utilizing the 
formalism of Majorana modes). We also present a DMRG-like algorithm to obtain the compressed 
correlation matrix directly from a hopping Hamiltonian. In addition, we discuss a slight variation of 
the procedure which leads to a simple construction of the multiscale entanglement renormalization 
ansatz (MERA) of a fermionic Gaussian state, and present a simple picture of orthogonal wavelet 
transforms in terms of the gate structure we present in this paper. As a simple demonstration 
we analyze the Su-Schrieffer-Heeger model (free fermions on a ID lattice with staggered hopping 
amplitudes). 


I. INTRODUCTION 

One of the streng ths of the density matrix renormaliza¬ 
tion group (DMRG)pJ21, and tensor network states in gen¬ 
eral, is that their power to simulate strongly correlated 
systems does not require the interactions to be weak. 
In fact, in fermion systems such as the Hubbard model, 
DMRG is more accurate for larger interactions. The ma¬ 
trix product state (MPS) representation of the wavefunc- 
tion, which DMRG implicitly uses, more efficiently com¬ 
presses the wavefunction when interactions are strong, 
due to lower entanglement in a real-space basis. 

In this paper, we introduce a new algorithm for ef¬ 
ficiently producing an MPS representation for ground 
states of noninteracting fermion systems. Why is this 
useful, when DMRG is most useful in the opposite 
regime? This would be a valuable tool in a number of 
situations. For example, a powerful and widely used 
class of variational wavefunctions for strongly interacting 
systems begin with a mean-field fermionic wavefunction, 
and then one applies a Gutzwilier projection to reduce or 
eliminate double occupancy.!^ It could be very useful to 
find the overlaps of a DMRG ground state with a variety 
of such Gutzwiller states to help understand and classify 
the ground state. Once one has the MPS representation 
of the mean field state, the Gutzwiller projection is very 
easy, fast, and exact, whereas in other approaches it usu¬ 
ally must be implemented with Monte Garlo. One might 
also begin a DMRG simulation with such a variational 
state, or in some cases with a mean field state with¬ 
out the Gutzwiller projection. Being able to represent 
fermion determinantal states as MPS’s in a very efficient 


way also opens the door to using DMRG ground states as 
minus-sign constraints in determinantal quantum Monte 
Carlo, in particular in Zhang’s constrained path Monte 
Carlo (CPMC) method^. In this case one would hope 
that, for systems too big for accurate DMRG, at least 
the qualitative structure of the ground state could be 
captured by DMRG, and then the results could be made 
quantitative with the Monte Carlo method. 

The basis of our approach shares ideas with DMRG. 
Matrix product state representations exploit a property 
of the state (low entanglement) to compress the informa¬ 
tion in the state. Fermionic Gaussian states (the gen¬ 
eral class of states which includes both fermion determi¬ 
nants, BCS states, and free fermion thermal states) are 
also compressible, as we will show. The properties of a 
Gaussian state are completely defined by its correlation 
matrix. For the case of a fermion determinant, the cor¬ 
relation matrix has eigenvalues which are either 0 or 1, 
i.e. they carry only a limited amount of information, in¬ 
dicating that the state can be compressed. In particular, 
one can perform an arbitrary single-particle change of ba¬ 
sis within the occupied states, or within the unoccupied 
states, without changing the determinantal state. Ten¬ 
sor network methods in the context of fermionic Gaus¬ 
sian states have been studied previously in the context 
of the multiscale entanglement renormalization ansatz 
(MERA)P and projected entangled pair states (PEPS)P, 
however here we present a simple and easily generaliz- 
able formalism and construction starting with an efficient 
method for forming the MPS of a fermionic Gaussian 
state. We also present a new and simpler method for 
obtaining a fermionic Gaussian MERA (GMERA), the 
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MERA of a fermionic Gaussian state, as a simple exten¬ 
sion. 

Our approach to producing the MPS of a fermionic 
Gaussian state also produces a compressed form of the 
correlation matrix itself, which we call a fermionic Gaus¬ 
sian MPS (GMPS), which might be useful in very differ¬ 
ent contexts where the single-particle matrices are very 
large. This compressed form expresses the N x N cor¬ 
relation matrix in terms of 0{BN) real angles which 
parametrize nearest neighbor rotation gates, where B <C 
N for states with low entanglement. The compressed 
form can be utilized directly. For example, ordinarily 
multiplying an arbitrary vector by the correlation ma¬ 
trix, which is not sparse, requires 0{N'^) operations, but 
by using the compressed form only 0{BN) operations are 
needed. For simplicity, the algorithm we introduce first 
utilizes the correlation matrix as the initial input. How¬ 
ever, in Appendix we present a DMRG algorithm in 
the single particle context, which we call fermionic Gaus¬ 
sian DMRG (GDMRG), that starts with a single par¬ 
ticle Hamiltonian matrix and outputs the ground state 
correlation matrix in compressed form as a GMPS at a 
greatly reduced cost compared to directly diagonalizing 
the Hamiltonian matrix, 0{B^N) as opposed to 0{N^). 
This algorithm exploits the close relationship between 
the correlation matrix and the density matrix of a many 
particle state, and many tensor network algorithms can 
similarly be translated into a single particle framework. 

The paper is organized as follows. Section [IT] gives a 
brief overview of fermionic Gaussian states and correla¬ 
tion matrices, including an introduction to the entangle¬ 
ment of these states. In Section [Till we give detailed de¬ 
scriptions of the new algorithms. Section HI A covers our 
new procedure for compressing a correlation matrix as a 
GMPS. Section |HI B| presents a variation of this method 
to obtain a GMERA. In Section [ill C| we give a brief in¬ 
troduction to how the GMERA gate structure relates to 
wavelet transforms. Section HID covers the procedure 
for turning the gates obtained from compressing the cor¬ 
relation matrix into a many-body MPS approximation of 
the Gaussian state. Finally, Section IV shows numerical 
results for the algorithms covered in the paper. 


II. BACKGROUND ON FERMIONIC 
GAUSSIAN STATES AND CORRELATION 
MATRICES 


Consider the Hamiltonian for a ID system of noninter¬ 
acting fermions 


N 

(1) 

* j=i 

where and a| are fermion operators and H = [Hij] 
is a symmetric matrix [H = H^). We assume that the 
Hamiltonian terms are local (so the matrix H is band- 
diagonal) . 


Diagonalizing the matrix H, we have H = UDW 
where U is orthogonal and D is diagonal such that 
Dkk' = ^k^kk'■ The Hamiltonian can then be put into 
diagonal form, 

N 

H = '^CkalcLk ( 2 ) 

k^l 

where the operators which create the single particle en¬ 
ergy eigenstates are 

N 

al = '^Uika\- ( 3 ) 

Assuming \i k < k' ^ the ground state is 

Nf 

IV'o) = it 4 1^) • (4) 

k=l 

where Np is the number of particles in the system. 

The correlation matrix is 

Np 

= = (5) 

fc=i 

The correlation matrix fully characterizes |^/'o) because 
all correlation functions, and therefore all observables, 
can be factorized into two-point correlators using Wick’s 
theorem. Note that the eigenstates of H are also the 
eigenstates of A (the same U that diagonalizes H also 
diagonalizes A). However, the eigenvalues of A are either 
1 (occupied) or 0 (unoccupied). The massive degeneracy 
of A means that we can make arbitrary changes of basis 
among the eigenstates of A as long as we do not mix 
occupied and unoccupied states. 

In our procedure, we will be interested in finding lo¬ 
calized eigenvectors of the correlation matrix which are 
(approximately) fully occupied or unoccupied. By ro¬ 
tating into the basis of these eigenvectors, we can lo¬ 
cally diagonalize the correlation matrix, which will lead 
to a compression of the state. These eigenvectors have 
eigenvalues near I or 0, which makes them (approximate) 
eigenvectors of the entire correlation matrix and therefore 
uncorrelated with the rest of the system. What makes it 
possible to find a localized eigenvector? 

The answer is the limited entanglement structure of 
the states we are interested in (ground states of local 
Hamiltonians). Gonsider the entanglement entropy of 
our fermionic Gaussian state, which can be calculated 
directly from the correlation matrix. Divide the system 
into an arbitrary subblock B oi B sites (with the cor¬ 
responding submatrix of A, which we call Ag) and the 
rest of the system. We would like to know how large of 
a block size B we need to find a localized eigenvector. 
If the matrix Ag has eigenvalues {n^} for b G B, with 
0 < rif, < 1, the entanglement entropy of the subblock B, 
Sb = — Tr[/5B log(/5B)] (where pg is the reduced density 
matrix of the state in subblock B), is 

SB{{nb}) = '^Si{nb) 
beB 


( 6 ) 
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(a) Eigenvalues and entanglement entropy. 



(b) Example eigenvectors. 


FIG. 1. Fig.j^a) shows the occupations Ub and corresponding 
entanglement Si{nb) from diagonalizing a block of B = 16 
sites in the middle of a system of free gapless fermions on 
N = 1000 sites at half filling. The minimum and maximum 
eigenvalues, ni and nie, differ from 0 and 1 by « 1.74x 10“^^. 
The eigenvalues closest to 1/2, 1/2 — ng = ng — 1/2 « 0.21, 
have entropies Si{ns) = Si{ng) ~ 0.60, which are close to 
the maximum of Si(l/2) = log(2) « 0.69. Fig. |^b) shows 
examples of eigenvectors from the same diagonalization. The 
eigenvectors with eigenvalues near 0 and 1, which contribute 
very little to the entanglement, are localized in the middle 
of the block, while the eigenvectors with eigenvalues closer 
to 1/2 which contribute most to the entanglement have large 
support on the edges of the block. 


where Si{nt) = -[n&log(nh) + (1 -nh) 1(^1 - n&)]. This 
expression has been shown elsewhere^nlil We show a 
simple, self contained derivation of it in Appendix [A| 
Note that Si{nh) vanishes for both n;, —>• 0 and Uf, —>• 1. 

The maximum amount of entanglement a block of size 
B can contain is when ni, = 1/2 for all b G B, so 
Sb < .Blog(2). This reflects a volume law entanglement 
in the “volume” B. However, ground states of ID local 


Hamiltonians have entanglement that is much smaller, 
either of order unity (if the system is gapped), or the 
entanglement grows as log(H) if the system is gapless. 
To avoid the volume entanglement, most of the block 
eigenvalues Ub must be exponentially close to 0 or 1. In 
other words, as soon as we make B big enough so that 
the entanglement begins to saturate, except for a possi¬ 
ble slow logarithmic growth, we should find at least one 
eigenvalue very close to 0 or 1. For gapless free fermions 
in ID on Af = 1000 sites, we show example eigenval¬ 
ues, eigenvectors, and corresponding entanglements of a 
block of i? = 16 in the middle of the correlation matrix 
in Fig. Even for gapless free Fermions, with a block 
size of only B = 16 we find many eigenvalues near 0 or 
1 (many localized eigenvectors). We use this observation 
next to develop algorithms to locally diagonalize correla¬ 
tion matrices and in the process find a very compressed 
form. 


III. ALGORITHMS 

A. Compressing a Correlation Matrix as a GMPS 


We begin the procedure by diagonalizing the upper 
left B X B subblock of a correlation matrix A of a pure 
fermionic Gaussian state. Assume that the state has 
some local entanglement structure, for example it is the 
ground state of a local Hamiltonian in ID. For now, we 
imagine our system has open boundary conditions. Let B 
be the group of sites 1,... H on the left end of the system, 
and Ag be the associated subblock of A. Also, let {ub} 
be the eigenvalues of Ag for 6 e H where 0 < Uf, < 1. 
(This constraint on the eigenvalues of the subblock fol¬ 
lows from the fact that both A and 1 — A are positive 
semi-definite.) 

We increase B until we find some Ub that is nearly 1 or 
0 within a specified tolerance, e.g. 10“®. The closer the 
eigenvalue is to 1 or 0, the more accurate the compres¬ 
sion, but a larger block size translates to more gates and a 
bigger bond dimension of the MPS we will form. In Fig.[^ 
we show the most occupied and unoccupied eigenvectors 
of Ag for H = 12 for a system of gapless free fermions in 
ID with N = 1024 sites. We see that H = 12 is sufficient 
to give deviations from occupancies of 0 or 1 to nearly 
machine double precision. The eigenvalues found in the 
bulk likely will not be as accurate, because states in the 
bulk will generally be more entangled than the ones on 
the edge. The smooth fall-off to zero at the right edge of 
the block is characteristic of these modes and is a result 
of diagonalizing the block on the left-most boundary of 
the system. The localized states we find here are least 
entangled with the rest of the system. This is in contrast 
to the dominant Schmidt states that are utilized within 
DMRG which have degrees of freedom that are localized 
at the edge of the block. 

The eigenvector v which is least entangled is also an 
approximate eigenvector of the total correlation matrix 
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FIG. 2. Examples of approximate occupied and unoccupied 
eigenvectors of A obtained from diagonalizing Ag where sub¬ 
block B are sites 1,... ,B. A is formed from the ground state 
oi H = + h.c.) for N = 1024 at half fill¬ 

ing [Nf = N/2). a block size of B = 12 is used. Eigen¬ 
vectors with highest (riocc) and lowest (nunocc) eigenvalues 
found from diagonalizing subblock B are shown. We find 
1 — riocc = 2.4 X 10“^® and riunocc = 7.3 x 10“^®, so the 
occupations are accurate to nearly machine double precision. 
1 — Wocc and riunocc should be equal at half filling (because 
of particle-hole symmetry), but are different in this case as a 
result of roundoff errors. 


A, i.e. AE « riiv. Any N x N unitary matrix that has v 
as its first column represents a change of basis that puts v 
on the first site. The associated transformation of A will 
make An = ni, and zero out the rest of row 1 and column 
1. The matrix of eigenvectors of Ag would produce such 
a matrix (expanding it to x Af by putting ones on the 
diagonal), but this B x B matrix does not translate well 
to many-particle gates to use in constructing an MPS. 

We now introduce gate/circuit diagrams which apply 
equally well to simple matrix manipulations of A and 
to many-particle tensor networks. The basic ingredient 
of the diagrams are two site nearest neighbor unitary 
gates. In Figure]^ we show the relation between a gate 
and a matrix. In the matrix interpretation, an “MPS” 
is just a vector with a dimension equal to the number of 
sites. In a later section we show how a gate is interpreted 
in the many-particle context of a tensor network. We 
consider nearest neighbor gates because these translate to 
fast MPS algorithms—typically, a non-nearest neighbor 
gate is implemented as a set of swap gates to bring the 
sites together, a nearest neighbor gate, followed by swaps 
to return to the original ordering of the sites, which is 
much slower than a single nearest neighbor gate. In the 
special case that the intermediate sites are in product 
states, i.e. bond dimension 1, non-local gates are also 
inexpensive, and we use these in our MERA algorithm. 

Returning to the task of moving the least entangled 
state V to the first site, a set of R — 1 two-site gates 
suffices. The first gate acts on sites (B — 1,B), and we 


V{03) 


/I \ 

1 

Vids) j 

1 

V V 


^(^3) 


/ cos 03 — sin 03 

sin 03 cos 03 


FIG. 3. Definition of a gate used throughout the paper. Ex¬ 
ample for A’ = 8 sites for a gate at site i = 4. Unless specified 
otherwise, circuits are in a direct sum space. We take the 
convention that multiplying a matrix from the top by a vec¬ 
tor corresponds to multiplying the matrix on the right by a 
column vector. 


label it Vb-i- In general, we take 


= ViOi) 


f cos 9i — sin 9i 
lsin0i cos0i 


( 7 ) 


We choose 9b-i = ia.n~^{vB/vB-i), where Vi is 
the component of the (un)occupied eigenvector 
of interest v. With this choice, Vb-i acting on 
= (ui ... vb -1 vb) sets the last component, 
vb, to zero, and produces a new value of vb-i —t 
v'b_i- In other words, we solve for 9b-i so that 
^Vb-1 = {vi ... VB-l vb)Vb-i = {vi ... v'g_^ 0). 
Next we rotate sites {B — 2, B — 1), with 9b-2 = 
tan“^(u^_j^/tiB- 2 ), and continue in this fashion. The 
action of all these gates on iF gives so they act to 
change the basis into the one containing v. 

We take Vjs = V{9b-i)V{9b-2) ■ ■ ■ This proce¬ 

dure is shown schematically for a simple case in Fig.j^a). 
We then apply the gates to A. The transformed corre¬ 
lation matrix VgAUg will have ni « 1 or 0 as the top 
left entry (and nearly 0 in the rest of the entries in the 
first row and first column). A schematic of this trans¬ 
formation is shown in Fig. |^b). We will call the first 
block Bi = B. We repeat this procedure for B 2 , sites 
2,...,B + 1, now simply ignoring the first site. For 
half-filled systems, the modes found are likely to alter¬ 
nate between occupied and unoccupied because occupied 
and unoccupied modes will generally be found in pairs 
when diagonalizing a block of the correlation matrix. Of 
course, B does not have to stay the same from one block 
to the next, and in general it is better to set it dynami¬ 
cally to make sufficiently close to 1 or 0. For the last 
blocks, B is decreased to the remaining number of sites. 
After N blocks, we will have approximately diagonalized 
A. 

The overall unitary transformation is U = 
VbiVb 2 ■ ■The matrix V decomposed into 
the 2x2 rotation gates {V{9i)} for N = 8 and R = 4 is 
shown in Fig. [^a). The N x N unitary approximately 
rotates our single particle basis from real space to what 
we refer to as the occupation basis, which is one of the 
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(1 0 0 0 ) 


v{e,) 



v{e2) 


{vi v'i 0 0) 

02 = tan-i (^) 


{Vi V2 V3 O) 


v{es) 


03 = tan 1 (^) 


= (t>l V2 V3 V 4 ) 


Av ~ niv, fii ~ 1 

(a) 











A' 









(b) 


FIG. 4. In Fig. j^a) we show schematically the procedure to obtain, given an approximate eigenvector v of the correlation 
matrix A, the set of local rotation gates that make up our compressed correlation matrix. The example shown is for a block 
size B = A and system size N = 8. Fig. |^b) shows that, by conjugating the correlation matrix by the gates obtained, the 
correlation matrix is approximately partially diagonalized. 


highly degenerate eigenbases of the correlation matrix. 
Conjugating A by C gives us a matrix V^AV that is 
nearly diagonal, with Np values on the diagonal close to 
1 corresponding to occupied modes and N — Np values 
on the diagonal close to 0 corresponding to unoccupied 
modes. In total, the procedure as described would 
require 0{BN) nearest neighbor rotations, where B is 
the largest block size needed for the desired accuracy of 
the representation of the correlation matrix. 

Writing the 2x2 rotations as gates is very convenient 
for understanding the matrix transformations, but more 
importantly it makes it easy to connect to many-body 
gates and to quantum circuits in general. As a quantum 
circuit, these gates have a slightly peculiar structure. Be¬ 
cause of how the diagonalizations overlapped, the circuit 
has a depth of 0{N). However, a vertical cut through 
the circuit only passes through 0{B) gates. This con¬ 
struction and gate structure is in a certain sense optimal 
if we limit ourselves strictly to circuits with local gates. 
If we want to represent a correlation matrix in a com¬ 
pact way with nearest neighbor gates, we would like to 
be able to represent arbitrary correlations in the system 


(correlations at all lengths), and in particular, correla¬ 
tions between the first and last site. In Fig. [gb), we 
show a circuit which cannot connect the first and last 
sites because its depth is less than N/2. Although our 
gate structure, shown in Fig. [^a), has a depth ~ N, in 
fact we can adjust our diagonalization procedure slightly 
to obtain a depth of ~ N/2 so our circuit can capture 
correlations of all lengths. This is done by beginning the 
diagonalization procedure from both the left and right 
side of the system until the blocks meet in the middle. 
This freedom in where to start the diagonalization is sim¬ 
ilar to the choice of gauge of an MPS. Choosing one gauge 
over another can be useful if we have already performed 
this procedure for a correlation matrix and want to per¬ 
form it again for another correlation matrix which is only 
locally different from the first one. If we choose the gauge 
center where the correlation matrix has changed, we only 
need to change a local set of gates. 

A generic local circuit of depth 0{N) contains 0{N'^) 
gates, and can represent an arbitrary N x N single¬ 
particle unitary change of basis. The low entanglement 
of physical ground states allows us to represent an N x N 
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FIG. 5. Fig.j^a) shows the overall gate structure obtained by the diagonalization procedure. These gates form the total N x N 
unitary V which approximately diagonalizes our correlation matrix A. By conjugating a diagonal matrix with the appropriate 
occupations of 0 or 1 found in the diagonalization procedure by this set of gets, we get an approximation for the correlation 
matrix. Fig.j^b) shows an example of the correlations allowed by representing the correlation matrix A with a diagonal matrix 
conjugated by a finite depth circuit of depth < N/2. The grey area (the “light cone”) represents sites where there can be 
nonzero correlations with the first site. For this circuit depth, there can’t be correlations with the last two sites. A circuit of 
depth > N/2 is required to allow for the possibility of nontrivial correlations across the entire system. 


matrix with 0{BN) one-parameter gates, with B N. 
For a gapless system, we know from conformal held the¬ 
ory that the entanglement of a subblock B of B sites 
varies as Sb ^ log(S). This means that we should be 
able to capture the entanglement of a critical system of 
N sites with a block size B ^ log(iV). If we hnd that 
B ~ log(iV), this means that our construction is roughly 
optimal. Fig. in Section |IV A| shows numerical evi¬ 
dence that this is indeed the case. 


B. Compressing a Correlation Matrix as a 
GMERA 

A MERA tensor networlJi^ can represent a ID criti¬ 
cal system using a constant bond dimension, unlike an 
MPS. In our MPS construction, this is reflected in that 
B ~ log(Af). However, we can adjust the diagonalization 
procedure slightly to obtain a MERA-like gate structure 
with a B which does not grow with N. The MERA for 
fermionic Gaussian states was hrst studied iiP, but was 
only used to study inhnite translationally invariant sys¬ 
tems and required a subtle optimization scheme. Here 
we will show a simpler construction only requiring the 
tools we have explained so far. 

We begin the procedure in the same way as we did 
for the GMPS, by diagonalizing the block corresponding 


to sites 1,..., B of the correlation matrix. Just as be¬ 
fore, for a large enough block size we find an occupied 
or unoccupied mode and rotate into the basis containing 
that mode with B — 1 local 2x2 gates. Next, instead 
of diagonalizing the block starting at site 2, we instead 
diagonalize the block corresponding to sites 3,..., H -I- 2, 
again finding an occupied or unoccupied mode and rotat¬ 
ing into that basis. The state at site 2 is “left behind”—it 
is not a low entangled state, so we cannot ignore it, but 
we leave it for a later stage of the algorithm. We con¬ 
tinue in this manner, diagonalizing blocks starting at odd 
sites of size B to obtain ~ BN/2 nearest neighbor gates. 
Approximately half of the modes are fully occupied or un¬ 
occupied and are projected out (meaning the associated 
rows and columns in the correlation matrix are ignored 
in later stages). The other half were left behind, and con¬ 
tinue as the sites of the next layer of the gate structure. 
By only trying to get N/2 unentangled modes in the first 
layer, the size of B does not need to grow with N, as we 
show below. 

The left-behind sites pass through to the next layer and 
are interpreted as a course-grained version of our original 
state on only N/2 sites. We repeat the same procedure 
for this new course-grained system of A/2 sites, starting 
by diagonalizing the subblock of the first 1,... ,B sites 
of the new course-grained lattice, finding an occupied or 
unoccupied mode of the course-grained system, and pro- 
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FIG. 6. An example of an alternative diagonalization scheme 
resulting in a MERA-like gate structure. Here we show a 
section of the first two renormalization steps, with 12 sites 
shown in the first layer and 6 renormalized sites shown in the 
second. A block size of i? = 4 is used. For this block size there 
are two layers of disentanglers and one layer of isometries per 
level of the MERA. 


jecting it out. Here, however, the gate we use to rotate 
into the basis of the (un)occupied mode are 2x2 nearest 
neighbor gates in the course-grained lattice, but are ac¬ 
tually next-nearest neighbor gates acting on the original 
lattice (if we project out every other site). Ordinarily, 
using next-nearest neighbor gates (or longer range gates 
at higher levels of the MERA) would be costly in the 
many-body case, requiring swap gates to make them ef¬ 
fectively nearest neighbor. However, the projected-out 
sites are now in product states, meaning that swapping 
does not require significant time. 

We repeat the above procedure of projecting out ev¬ 
ery other effective site and course graining to a lattice 
of half the size. All of the sites will be projected out 
after this course-graining is repeated 0 (log 2 (fV)) times. 
Fig. [6] shows the first two layers of the resulting gate 
structure, which looks like a MERA with B — 2 layers of 
nearest neighbor 2-site disentanglers and a layer of near¬ 
est neighbor 2-site isometries. The total number of gates 
in the construction is ~ B{N/2 -|- 7V/4 -|- ... -I- 1) = BN, 
the same gate count for a fixed block size B as for the 
GMPS. We call this gate structure, which like our GMPS 
construction is a compression of an x correlation 
matrix into ^ BN gates, a fermionic Gaussian MERA or 
GMERA. In this figure, triangles denote projections onto 
the appropriate occupation found. These are the modes 
that are uncorrelated with the rest of the sites and can 
be ignored in the next layer. Some extra gates will be 
required to project out the leftover sites at the right end 
of the system (not shown in Fig. [^, and there is some 
flexibility in how to do this which will change the accu¬ 
racy of the compression slightly. For example, one could 


use a gate structure similar to the GMPS construction 
to project out all of the leftover sites at the end. 

How does the block size B of the GMERA compare to 
that in our GMPS construction? We show numerically 
in Section IV B| that for a simple gapless Hamiltonian 
the GMERA does indeed produce accurate results with 
a block size B = 0(1), independent of the system size, 
making it much more efficient in the large N limit. 


C. Discrete Wavelet Transforms and Fermionic 
Gaussian MERA 

We would like to point out the similarity between the 
MERA gate structure and orthogonal wavelet transforms 
(WT), such as the WTs that produce the well-known 
Daubechies wavelet^ ^^^t l Qf course, the development of 
wavelets has not been in a many particle context, and, 
for now, we restrict ourselves to the matrix interpretation 
of the diagrams. For compact wavelets, an orthogonal 
wavelet transform is a local unitary transformation. It 
is not usually represented in terms of two-site gates, but 
this representation turns out be be particularly conve¬ 
nient. To be specific, we start with the simplest nontriv¬ 
ial WT, the D4 Daubechies WT. This WT is defined by 
four coefficients {oj} for j = 1,..., 4 which characterize 
how the D4 scaling function is related to itself at differ¬ 
ent scales through the equation s{x) = ajV^s{2x—j). 
The matrix form of the WT is given by 


/ai 

02 

03 

04 

0 

0 

0 

(I4 

-03 

02 

-Ol 

0 

0 

0 

0 

0 

Ol 

02 

03 

(14 

0 

0 

0 

04 

-03 

02 

-Ol 
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V •■•/ 

The {uj} are carefully chosen to ensure orthogonality be¬ 
tween scaling functions centered at different sites, and to 
make the scaling functions have desirable completeness 
properties. For example, linear combinations of the D4 
scaling functions centered at different sites, {s{x — k)} for 
integer k, fit any linear function, so the resulting coeffi¬ 
cients are = (1 -I- -x/S, 3 -I- -x/S, 3 — -x/S, I — -x/S)/(4-\/2). 
The orthogonality requirement results in nonlinear equa¬ 
tions to solve for the {aj} which becomes complicated 
for higher order. The second row of the matrix gives 
the coefficients that produces wavelets, designed to rep¬ 
resent high momentum degrees of freedom. In terms of 
our MERA procedure, the wavelets are left behind, while 
the scaling functions propagate to the next level. 

The D4 WT has a very simple gate structure, identi¬ 
cal to our MERA structure with B = i, shown for two 
layers in Fig. In each horizontal layer of gates, all 
gates have the same angle. The D4 WT is specified by 
only two angles, 6 i for the bottom layer and O 2 for the 
next. Higher order WTs of this type (e.g. D6, D8, etc.) 
correspond to larger B. (For example, the D6 WT looks 
like Fig. 6). Given the angles, one gets the {aj} by set¬ 
ting all the top values of the circuit to zero except a I 
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FIG. 7. Here we show an example of a discrete wavelet 
transform written in the gate notation introduced in this 
paper. We show the D4 wavelet, which corresponds to a 
fermionic Gaussian MERA with one layer of disentanglers 
and one layer of isometries per layer. wi and Si (w 2 and 
S2) label wavelet and scaling functions for the first (sec¬ 
ond) layer. Taking 61 = tt/G and 62 = 57r/12, we repro¬ 
duce the conventional scaling coefficients for the D4 WT, 
= (fli, 02,03,04) = (1 -f %/3, 3-1- %/3, 3 — %/3,1 — %/3)/(4%/2). 


on one site and applying the 2 x 2 rotations in the layers 
below. The support of the scaling functions is made obvi¬ 
ous using the gate structure, as there will be 2 L nonzero 
values at the bottom of the circuit for L layers of gates. 
For the D4 WT, one finds that 9i = tt/G and 62 = 5tt/12 
reproduces the D4 {aj}, up to a trivial reversal of the co¬ 
efficients. (A single layer with 6*1 = 7 r /4 gives the trivial 
Haar wavelets, which have been used previously as a ba¬ 
sis for transforming fermionic Gaussian states by qP.) 
The scaling functions at the larger scales are found by 
performing the same transformation of L layers of gates 
on the scaling functions found at the previous scale. 

In Fig. I^we show how scaling coefficients {a^} come 
from the gate structure, applying a vector to the top of 
the circuit with 1 at the site of a scaling function and 
O’s elsewhere. In simple wavelet treatments, the wavelet 
coefficients are obtained from the scaling coefficients {oj} 
as {{—iy~^a 2 L-j+i} for j = 1,...,2L. Here, they are 
obtained by shifting the location of the 1 at the top of 
the circuit, but we can show in general that this gives the 
same result. This is done by noting that the shift of the 
1 to get the wavelet coefficients looks like a swap at the 
top of the circuit. We can “pull through” this swap by 
conjugating each layer of the WT with a transformation 
that reverses the order of the sites. This conjugation also 
negates the angles in the circuit. It leaves a site reversal 
at the bottom of the circuit, reversing the order of the 
coefficients. The angle negation negates the sine terms, 
leading to the same coefficients except with every other 
one negated, since every other site will have an even or 
odd number of sin( 0 i) multiplied together. 

Given an arbitrary set of {oj}, we can use the same 
procedure that brought v to the first site in our GMPS 
procedure to find all the angles defining the WT, i.e 
V = a. Thus, any compact orthogonal WT of this gen¬ 



fa) Scaling coefficients from gate structure. 
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and Ci = cos{9i), Si = sin(0i) 


(b) Gate structure in (a) written in terms of matrices and vectors. 



(c) Wavelet coefficients from gate structure. 



FIG. 8. Here we show explicitly how to obtain the scaling 
and wavelet coefficients of the D4 WT from the circuit con¬ 
struction. Taking 9i = tt/G and 62 = jVl, in (a) and (b) we 
reproduce the conventional scaling coefficients for the D4 WT, 

aF = (tti, 02 , as, ay = (1-|- \/3, 3-1- \/3, 3 — \/3,1 — \/3)/(4\/2), 

up to a trivial reversal in the order. In (c) with the same 
choice of angles we reproduce the conventional wavelet coef¬ 
ficients (a 4 , —as, 02 , —ai), again up to a trivial reversal and 
sign. 


eral type can be represent by a simple gate structure. 
Because wavelets are much easier to understand than 
generic many particle wavefunctions, the connection be¬ 
tween MERA and wavelets may help provide intuition 
that helps one understand MERA. 


D. Forming the Many-Body MPS from the GMPS 
(or GMERA) 


For a number conserving real Hamiltonian H, 
the many particle unitary gate Vi corresponding 
to the single particle rotation V), in the basis 

{|fl) , a| \n) , aj+i |H), ajaj+i |fl)}, is 


= [V{9^)] 


/I 0 0 0\ 

I 0 cos 9i sin 9i 0 1 
I 0 — sin 9i cos 9i 0 I 

Vo 0 01 / 


(9) 
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FIG. 9. Tensor diagram showing the structure of gates {t^;} for i = 1,..., N — 1 obtained in our procedure and how they 
contract to form an MPS. Here we show a system with N = 8 sites and a block size 5 = 4. The diagram on the right 
shows that once the sites are rotated into a basis where one of the modes is occupied or unoccupied (generally with some 
alternating pattern), the fully occupied or unoccupied modes can be projected out. The transformations {Vb^}, including the 
projections, can be directly interpreted as the tensors composing the MPS representation of our many-body state if we do an 
exact contraction, or we can apply them as a set of gates as explained in the text. 


This reinterpretation of the gates is the only change need 
to make our matrix gate structures act on the many par¬ 
ticle Hilbert space. 

Say we have compressed the correlation matrix of a 
pure fermionic Gaussian state as a GMPS. To create the 
MPS representation of this state, we begin with a prod¬ 
uct state, with each site being occupied or unoccupied, 
with the occupations given by n/^ obtained in our diag- 
onalization procedure (set to 1 or 0 for « 1 or 0). 
We then apply, one by one, all of the nearest neighbor 
gates {Vi} (the many-body gates corresponding to the 
gates {Vi} obtained with Eq. in the opposite order in 
which they were obtained with our diagonalization proce¬ 
dure. The repeated application of gates is similar to the 
time-evolving block decimation (TEBD) algorithrrfi^ or 
the time dependent DMRG algorithrrP^^ but the pattern 
of gates and ordering is different. We apply the two body 
gates by moving the center of the MPS to the location of 
the gate, contract the gate with the two relevant tensors 
in the MPS, and then form the new MPS by perform¬ 
ing a singular value decomposition (SVD), with possible 
truncation of states by throwing out states with small 
singular values. 

We can also form the MPS from our GMERA construc¬ 
tion in a similar manner. However, instead of starting 
with a full product state, we start with the gates at the 
top of the MERA and work our way down, including only 
the sites that have been touched by a gate at that level 
or above. When a site is added, it starts as a completely 
occupied or unoccupied state, and is immediately mixed 
with another site by a gate. The number of sites involved 


roughly doubles with each layer, and after 0 (log 2 (iV)) 
layers of the MERA we have our MPS approximation for 
the entire system. Again, we can truncate as needed by 
throwing out low weight states after the SVD as we work 
our way down. 

Returning to the MPS construction, the tensors of the 
MPS could also be constructed directly by contractions 
of the gates as shown in Pig.[^ In this diagram the small 
black and white triangles signify projectors onto the ap¬ 
propriate occupations found, while the thick lines signify 
combined internal indices which form the internal bonds 
of the MPS. Prom this perspective it is easy to see that 
picking a block size B for diagonalizing the correlation 
matrix would correspond to an MPS with a bond dimen¬ 
sion of X = 2^~^. We find it simpler and more efficient 
to apply the gates layer by layer instead of this method. 
Layer by layer, it is natural to truncate the MPS with 
SVDs during the construction, and this can lead to an 
MPS with a smaller bond dimension than for the 

required accuracy. The SVD truncation takes one out 
of the manifold of Gaussian states, where the greater 
freedom for a fixed bond dimension allows one to find a 
state which is closer to the desired Gaussian state than 
one could within the Gaussian manifold. However, one 
should pick a block size so that 2 ^~^ is as close to the 
target bond dimension as possible. 

We can adapt our circuits to complex quadratic Hamil¬ 
tonians, where the gates are of the same form but the 
2 x2 submatrix rotating the singly occupied subspace is 
a general matrix in SU(2) parameterized by two angles. 
Even more generally, we can extend this procedure to 
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quadratic Hamiltonians with pairing terms to compress 
BCS states, where the gates required are not just num¬ 
ber conserving but general parity conserving gates (so 
they involve mixing of unoccupied and doubly occupied 
subspaces of the 2 sites of interest). This matrix would 
in general be parameterized by 5 angles (one matrix in 
SU(2) rotating the singly occupied subspace, one matrix 
in SU(2) rotating the empty and doubly occupied sub¬ 
spaces, and a relative phase). This form of gates has 
been studied previously in the context of classically sim¬ 
ulating quantu m circ uits using the matchgate formalism; 
see for exampld ^^bs i 


IV. NUMERICAL RESULTS 


cq 
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O 
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Energy gap 


Here we show numerical results for the algorithms we 
presented. In order to study systems that are both 
gapless and gapped, we study a simple model, the Su- 
Schrieffer-Heeger modeP^. This is a model of ID spin¬ 
less fermions hopping on a lattice with staggered hopping 
amplitudes, ti and t 2 - The Hamiltonian is 

iV-1 

2 

Hssh = ^ (tiaL-i«2* + t2a\^a2^+l + h.c.). ( 10 ) 

We will take ti = —t (l -I- |) and <2 = (l ~ f)- The 
model has an energy gap in the bulk between the ground 
state and first excited state of A = 2|5|t in the thermody¬ 
namic limit {N —> 00 ). With open boundary conditions, 
it can contain exponentially decaying zero energy modes 
localized on the ends of the chain. 


A. Results for Compressing a Correlation Matrix 
as a GMPS 


FIG. 10. Block size required to obtain a relative error in the 
total energy of less than 10“® as a function of the calculated 


energy gap (in 

f) 1 « 

units of t) for N = 128 sites. 


I ^ ^ 1 ^ 1 ^ 1 


• 

0.01 


bJD 

• 

g 0.001 


.3 

• 

g 0.0001 

• ^ 

b 

: : 

^ le-05 

• g 


: : 

(g le-06 

• 


1 = 

le-07 



_^^_I_^_I_L 

2 4 6 8 


Block size B 


We start with a simple test of obtaining the GMPS 
compression of the ground state correlation matrix of the 
SSH model for TV = 128 lattice sites for various energy 
gaps at half filling {Np = N/2). We analyze the range 
of 6 from 0 to 2. The ground state for 5 = 0 is (approx¬ 
imately) gapless while for i5 = 2 it is fully gapped (the 
chain uncouples). Fig. 10 shows the block size required 
to obtain a GMPS with a relative error in the total en¬ 
ergy of less than 10“® as a function of the calculated en¬ 
ergy gap. The exact ground state energy and energy gap 
are calculated by diagonalizing the hopping Hamiltonian 
Hssh- This corresponds to the accuracy of the MPS 
representation of the ground state if the GMPS written 
with many-body gates is contracted with no further trun¬ 
cation of the MPS, so a GMPS block size B corresponds 
to an MPS of bond dimension y = (which is why 

the block size remains constant for intermediate energy 
gaps). The plot shows, as expected, that the block size 
required decreases as the energy gap is increased. Fig.[TT] 
shows, for the case 5 = 0 (where the energy gap, due to 


FIG. 11. Relative error in the total energy as a function of 
the block size B for N = 128 sites and 5 = 0. 


the finite size, is 0.146088t), the relative error in the en¬ 
ergy as a function of the block size. 


Fig. 12 shows examples of the modes obtained with 
the procedure, both filled and unfilled, for a small gap 
and a larger gap. The modes are seen to be localized 
for the case of the larger gap, and extend throughout the 
system for the smaller gap. The unfilled modes follow the 
same decay as the filled modes but oscillate more, since 
they are above the Fermi sea and are therefore higher in 
energy. Fig.j^shows for the same two gaps the deviation 
in the eigenvalues Uk from 0 or 1 obtained during the 
diagonalization procedure. For the case of the larger gap, 
this error saturates to its maximum quickly for modes 
near the middle of the system, while for the smaller gap, 
the error increases more slowly due to the longer range 
correlations. 

In Fig.l^we analyze the block size scaling with system 
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Site 

(a) 5 = 0.4 (energy gap 0.8061354) 



FIG. 12. Examples of occupied and unoccupied modes 
found in the diagonalization process. Fig. |12[ a) shows occu¬ 
pied/unoccupied modes for <5 = 0.4 (energy gap « 0.806135t). 
Fig. [^b) shows occupied/unoccupied modes for 5 = 0 (en¬ 
ergy gap Rs 0.146088). 




Mode k 

(b) <5 = 0 (energy gap Ri 0.1460884) 

FIG. 13. Examples of deviations in occupations at the end of 
the diagonalization procedure for N = 128 sites. Fig. Ela) 
shows errors in the occupations for S = 0.4 (energy gap « 
0.8061354). Fig. |13[ b) shows errors in the occupations for 
S — 0.0 (energy gap « 0.1460884). 


size N for the gapless case {6 = 0). As we expect from 
arguments about entanglement made at the end of Sec¬ 
tion III A the scaling is found to be S ~ log(A4). This 
is the expected scaling for a critical ID system. We can 
see that with this procedure we can analyze very large 
systems, up to IV = 2^® = 65536 sites, even for gapless 
free fermions. To avoid storing correlation matrices this 
large, we begin with a very accurate compressed corre¬ 
lation matrix as a GMPS using the GDMRG algorithm 
presented in Appendix [B| With GDMRG, we begin with 
a state with a relative error in the energy of < 10“^®. 
For N = 65536 this requires a block size of i? = 22. We 
then obtain the local correlation matrix for the block we 
are interested in using the gates from this accurate com¬ 
pression, and use it to obtain a less accurate compression 
with a smaller block size. This procedure should lead to, 
for a given block size, a more accurate overall state than 
one that would be obtained directly from GDMRG, be¬ 


cause GDMRG optimizes the energy which only depends 
on very local correlations. 


B. GMERA Results 

Here we present results for compressing a correlation 
matrix as a GMERA using the procedure presented in 
Section |IIIB| We show the relative error in the energy 
for increasing number of sites for H = 10 in Fig. We 
see that for this block size, the error stays below 10“® 
for systems up to fV = 2^^ = 16384 and in fact appears 
to saturate at high number of sites (the change in the 
relative error in the energy approaches 0 for larger sys¬ 
tem sizes). This is in stark contrast to the GMPS, where 
a block size B ^ log(iV) was required to obtain a fixed 
accuracy, as shown in Fig. Instead, the GMERA ob¬ 
tains the same accuracy with constant block size B as 
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FIG. 14. Block size B needed for a relative error in the energy 
of < 10~® as a function of number of sites N for spinless, gap¬ 
less fermions with open boundary conditions at half filling. As 
expected from arguments about the entanglement of a critical 
system, we find B ~ log(A^), tested up to = 2^^ = 65536 
sites (note the log scale on the x axis). To study systems 
of this size and avoid the 0{N^) diagonalization of the hop¬ 
ping Hamiltonian, we obtain the correlation matrix using the 
GDMRG algorithm as explained in Appendix [B| 
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FIG. 15. Relative error in the energy for the proposed 
GMERA construction for increasing number of sites for a 
block size B = 10. The system analyzed is the ground state 
of free fermions hopping on a lattice with open boundary con¬ 
ditions. All errors are below 10“®. As expected for a MERA, 
the error is seen to saturate for large N, indicating a fixed 
block size is sufficient to obtain an accuracy < 10“® up to 
very large system sizes. 



shown in Fig. The GMPS obtains the given accuracy 
with the same or smaller block size up to iV = 512, after 
which it requires a larger block size than the GMERA to 
obtain the same level of accuracy. As we mentioned ear¬ 
lier, this is made possible partially because the GMERA 
structure involves nonlocal gates. 


FIG. 16. A plot of the time to form the MPS approximation 
of gapped and gapless free fermion ground states at half fill¬ 
ing as a function of sites N using gates from a GMPS. The 
bond dimensions are chosen large enough such that the rel¬ 
ative errors in the energy of the MPS are below 10“®. The 
block size of the GMPS used to form the MPS are the mini¬ 
mum required to obtain a GMPS with a relative error in the 
energy of 10~®. A cutoff in the singular values of the SVD 
of 10“^^ was used when applying the gates to form the MPS 
using the method described in Section |mg For the gapped 
case, the SSH model with <5 = 0.1 is used, corresponding to 
an energy gap of A ~ 0.2t (exact as A —>■ oo). 


C. GMPS to Many-Body MPS Results 


Plots of the time it takes to form the MPS of the 
ground state of a gapless free fermion system for up 
to A = 1024 sites using the method presented in Sec¬ 
tion |III D| are shown in Fig. As expected, the time 
it takes for a gapless system is polynomial in the sys¬ 
tem size A, while it is approximately linear in A for a 
gapped system. The SSH model is used with S — 0.1 or 
an energy gap A « 0.2t. The time to form the gapless 
ground state is only a modest polynomial in A, ^ A^ °^, 
while as we expect from arguments about entanglement 
the time to form the gapped ground state is very nearly 
linear in A, ^ A^ °^, because the block size and bond 
dimension required to obtain the specified accuracy are 
constant for all A shown {B = 8 and x = 55). With this 
method, a gapless ground state of A = 1025 sites with a 
relative error in the energy of < 10“® can be formed on 
a laptop in only ~ 90 seconds. 

An interesting point to emphasize is the quality of the 
compression. The GMPS for the gapless ground state on 
A = 1025 sites requires a block size of H = 11 to obtain 
a relative error in the energy of < 10“®. Naively, turn¬ 
ing these gates into many-body gates and contracting 
the network (forming the MPS directly from the GMPS 
with no truncation) as explained earlier leads to a bond 
dimension of the MPS of y = = 2^° = 1024. How¬ 

ever, applying the gates as described and using a cutoff 
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of the singular values of 10“^^ leads to the formation of 
an MPS approximation of the fermionic Gaussian ground 
state, still with a relative error in the energy of < 10“®, 
with a bond dimension of only x = 364. This is a result 
of the fact that our GMPS only explores the manifold of 
fermionic Gaussian states limited to the specified block 
size. On the other hand, the MPS approximation of the 
Gaussian state we form from this GMPS is able to ex¬ 
plore the entire manifold of MPS’s up to the allowed bond 
dimension (and particle number if symmetric tensors are 
used, as we do here), so through the SVD we are able 
to compress the state quite efhciently beyond what we 
initially might expect. 


V. CONCLUSION 


We have presented an efhcient, numerically stable, and 
controllably accurate way to compress a correlation ma¬ 
trix into a set of 2 x 2 unitary gates. From these gates, 
we have also presented a method for easily and efficiently 
forming the MPS approximation of a fermionic Gaus¬ 
sian state. We explained the procedure in detail for the 
ground state of a generic number conserving Hamilto¬ 
nian. We then presented results for the SSH model, a ID 
chain of fermions with staggered hopping. We showed 
examples of the accuracy and block sizes needed for dif¬ 
ferent gaps of the model. We hope this method can be 
used as a simple, efficient and reliable procedure for di¬ 
rectly preparing many states of interest, either by creat¬ 
ing starting states to aid DMRG calculations or prepar¬ 
ing a particular ansatz as an MPS. We also presented one 
example of how the procedure can be modified to obtain 
different gate structures, in this case one that is related 
to the MERA. However, there are other possibilities to 
be explored, such as gate structures more directly suited 
for systems with 2 spatial dimensions, periodic boundary 
conditions, as well as how the method might be applied 
to study thermal fermionic Gaussian states. In addition, 
we presented how discrete wavelet transforms can be de¬ 
scribed very simply with the gate structure notation we 
introduced in this paper. 

The method is easily generalized to cases beyond the 
one presented here. As we touched upon earlier, it can be 
generalized to the case of BCS states, the ground states of 
hopping Hamiltonians that include pairing terms. In this 
case, the correlation matrix in the Majorana basis can be 
written in terms of an anti-symmetric matrix which can 
be approximately block diagonalized by ~ 5BN local 2 x 
2 rotation gates, which are turned into nearest neighbor 
parity-conserving many-body gates. The case of spinless 
fermions was presented, but spinful fermions are a simple 
generalization. 
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Appendix A: Calculation of the Entanglement 
Entropy of a Fermionic Gaussian State 

In this section we give a simple, self-contained deriva¬ 
tion for Eq. the entanglement entropy for a block of 
a free fermion system. Assume the block of interest B is 
the first B sites. 

Gaussian states have expectation values that obey 
Wick’s theorem. This means that the expectation value 
of any operator contained within the block is specified if 
we know subblock B of the correlation matrix, Ag. This 
implies that the many-body density matrix of the block 
Pb is also uniquely specified by Ag. The entanglement en¬ 
tropy on block B, defined as Sb[pb\ = — Tr[/3g log(/3g)], 
does not change under general unitary transformations 
within the block. Thus, we can perform the single parti¬ 
cle unitary transformation of basis that makes Ag diag¬ 
onal, with diagonal entries rn, = (^a\a}}j for 6 S 1,..., H. 
The rib uniquely specify the reduced density matrix of 
the block, so the entanglement is a universal function of 
{nb}: 

Sb = SB{ni,...,nB). (Al) 

In fact, the details of the system outside the block are 
irrelevant. For example, different systems with different 
numbers of sites outside the block can have the same Sb 
as long as their {rib} are identical and the system is a 
Gaussian state. Thus to evaluate the function 5'g, we 
can choose a simple system in which to evaluate it rather 
than using the actual system of interest. 

Let’s first consider a block with only one site (H = 1). 
We would like to know the universal function We 

choose a two site system containing a single particle, with 
normalized wavefunction 

IV’) = -I- VI - riial) |H). (A2) 

The correlation matrix is 

f ni ^/nl{l-nl)\ , . 

W^i(l - ^i) 1 - ui ) 
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which has the required block correlation matrix Ai = 
(ni). The Schmidt decomposition of I'i/') is 

|V>) = ^/r^(|0))(aj |0)) + |0))(|0)) 

= VT^|0)|1) + Vrrr|l)|0) 

where |1) is the occupied state of the corresponding site. 
From this we see that the reduced density matrix for site 
1, Pi = Tt 2[\^) i'lpD, is 

Pi = (1-m) |0) (0|+ni |1) (1| 

(A5) 

= (1 - ni)(/i - ni) + nini 

and 

-S'i(ni) = - Tr[/5i log(/5i)] 

= -(1 - ni) log(l - ni) - m log(ni). ^ 

For the system B with block size B > 1, we can choose 
the system to be of size 2B and for each site in the block 
associate one site outside the block. The Gaussian state 
is the product state of the single particle states living 
on a pair, each identical in form to the B = 1 state, 
with B total particles. This system has no correlations 
or entanglement between these pairs. This means that 
the entanglement is the sum of the entanglement of each 
pair. Thus 

SBi{nb}) = '^Si{nb) (A7) 

beB 

which is identical to Eq. Note also that the overall 
reduced density matrix of the block is the product of the 
single site density matrices given in Eq. |A5| 

An alternative argument can be made to derive the 
same equation which avoids the introduction of a con¬ 
trived environment. We could have taken as an ansatz 
that the reduced density matrix on B, pg, is the product 
of the single site reduced density matrix we derived in 
Eq. |A5[ in other words pg = ®b^BPb- We can show that 
this is in fact the unique reduced density of the state we 
are interested in if we can show that it reproduces the 
correct correlation matrix of our state and is a fermionic 
Gaussian state (that it obeys Wick’s theorem). Both of 
these are easy to show explicitly. Once we verify that this 
is indeed the correct reduced density matrix of our state, 
we can calculate the entanglement entropy directly with 
Sb = - Tr[pg log(pg)] = -EhGBTr[pt,log(p6)], which 
matches Eq. 


Appendix B: GDMRG, an Algorithm to Obtain a 
Compressed Ground State Correlation Matrix as a 
GMPS 

Here we describe fermionic Gaussian DMRG 
(GDMRG), a DMRG-like algorithm in the single 
particle context. The algorithm is an efficient method to 
directly obtain all the angles specifying the compressed 
correlation matrix as a GMPS without ever needing to 



(b) 


FIG. 17. Fig. |17[ a) shows an example of an effective Hamil¬ 
tonian centered at site 4 for the GDMRG algorithm. The ex¬ 
ample is for N = 8 sites and a block size of R = 3. Here the 
center is only one site, but could be more to improve conver¬ 
gence just like in the standard DMRG algorithm. Fig. 
shows, for a sweep to the right, how to obtain the new left 
block from the previous effective Hamiltonian. 


express the matrix in uncompressed form. The ground 
state GMPS of a hopping Hamiltonian on N sites is 
calculated with a cost of only 0{B^N), where B is the 
block size of the GMPS (which determines the accuracy 
of the compression and depends on the entanglement of 
the ground state). 

Imagine that we start with a hopping Hamiltonian H 
on a lattice of N sites, and we would like to obtain the 
GMPS with a block size B that minimizes the energy of 
H. We begin with a random starting GMPS. Just like in 
the DMRG algorithm, we form an effective Hamiltonian 
centered at a site with a left and right block, which we 
show in Fig. Oa). Say that we start on the left side 
of the lattice and begin sweeping right. Our GMPS will 
start gauged to the left. For a single-site GDMRG, our 
center block is only one site, but we could use a larger 
center to decrease the number of sweeps required for con¬ 
vergence. In practice for free gapless fermions we find 
that a single center site works quite well. The first step 
is to diagonalize the 2B — 1 site effective Hamiltonian, 
and obtain the effective correlation matrix Aeff. Using 
this Aeff, we diagonalize the first B x B block and, for 
a large enough H, find a fully occupied or unoccupied 
mode. Just as described in Fig. we find the B — 1 
nearest neighbor 2x2 gates that rotate Aeff into the ba¬ 
sis containing this mode, partially diagonalizing it. These 
gates form the first block of the GMPS. 
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Next we would like to move the center to the right so 
that we can obtain the next block of the GMPS. Because 
the compression is a unitary transformation, we can start 
moving the center to the right by undoing the gates in 
the block of the GMPS to the right of the center. This is 
step is in contrast to ordinary DMRG where a sequence 
of right blocks are stored and are called from memory 
when needed. We then obtain the effective Hamiltonian 
for the next site using the block of the GMPS we obtained 
from the previous effective Hamiltonian of the sweep, as 


shown in Fig. ©b). We repeat this process until we 
reach the end of the lattice, completing our first sweep. 
We continue sweeping back and forth until the energy is 
sufficiently converged. We use this algorithm to obtain a 
very accurate correlation matrix for systems up to = 
2^® = 65536, from which we obtain the GMPS in Fig. 

For N = 65536 sites to obtain a correlation matrix with a 
relative error in the energy of less than 10“^*^, we require 
a block size of H = 22 and 14 sweeps (where a single 
sweep is from left to right or right to left). 
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